Problem statement: In this project, we tried to make classifications of the trails based on their length and change in elevation. We also made one function that would simplify the process of classification. To generalize this analysis, we used boundary, buildings, contours_3m, and trails layers from macleish package and contours250k layer from MassGIS.

#devtools::install_github("beanumber/macleish")
library(tidyverse)
library(sf)
library(macleish)
library(leaflet)
library(dplyr)
library(lwgeom)
library(ggplot2)
elevation <- st_read("contours250k/CONTOURS250K_ARC.shp")
## Reading layer `CONTOURS250K_ARC' from data source `/Users/ChristinaLyu/Desktop/SDS192/miniProject3/contours250k/CONTOURS250K_ARC.shp' using driver `ESRI Shapefile'
## Simple feature collection with 50395 features and 3 fields
## geometry type:  LINESTRING
## dimension:      XY
## bbox:           xmin: 33874.42 ymin: 778055.9 xmax: 325713.8 ymax: 959744.8
## epsg (SRID):    NA
## proj4string:    +proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000 +datum=NAD83 +units=m +no_defs
plot(st_geometry(elevation))

# add data source

We first plotted the elevation map of Massachusetts illustrates the number of feet or meters the state rises above sea level. Since the graph plots contours on the map, the denser the contour lines, the higher the place above the sea level. As you can see, the highest point in Massachusetts is in the northwestern part of the state, in Berkshire County, while the eastern third of Massachusetts is relatively lower than the other area.

elevation <- st_transform(elevation, 4326)
elevInter <- 
  elevation %>%
  st_intersection(macleish_layers[["boundary"]])

ggplot() + 
  geom_sf(data = elevInter, aes(color = factor(CONTOUR_FT), colors = "Set1")) + 
  geom_sf(data = macleish_layers[["trails"]], aes(color = factor(name), alpha = 0.8, lwd = 1.2)) +
  theme(legend.text=element_text(size=15))

leaflet() %>%
  addTiles() %>%
  addPolylines(data = st_geometry(elevInter), weight = 2) %>%
  addPolylines(data = st_geometry(macleish_layers[["trails"]]), color = macleish_layers[["trails"]]$color)

To rate the difficulty of the trail, we want to see the change of elevation. We made a contour map of the MacLeish Field Station, using different colors of lines to represent different contours and plotted the trails on the contour map. In this graph, each intersection of contour lines and trails represents a 30 ft change in elevation of that trail. To see intersections more clearly, we changed the alpha and the size of the lines so that pathways can stand out from the contour lines. To make an interactive visualization of elevation and trails, we created a leaflet map as well.

st_transform(elevation, 4326)
elevation %>%
  st_intersection(macleish_layers[["trails"]])
#change the color name in Macleish to lower case so that r can understand it
newTrails <- macleish_layers[["trails"]] %>%
  mutate(colorName = tolower(color))

markers <- 
  macleish_layers[["trails"]] %>%
  group_by(name) %>%
  filter(geometry == min(geometry))
#interactive map with buildings and the trails in Macleish
m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = macleish_layers[["buildings"]], popup = ~name) %>%
  addPolylines(data = newTrails, weight = 3, color = newTrails$colorName) 
m
leaflet() %>%
  addTiles() %>%
  addPolygons(data = macleish_layers[["trails"]], label = ~name) %>%
  addPolylines(data = newTrails, weight = 3, color = newTrails$colorName) %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Topo") %>% 
  # ...and a satellite imagery map layer
  addProviderTiles("Esri.WorldImagery", group = "Sat") %>%
  # Next you want to add the control to switch between the two
  addLayersControl(baseGroups = c("OpenSreetMap", "Topo","Sat"),
                   options = layersControlOptions(collapsed = F))

To rate the difficulty of the bike trails, we also want to look at the length of each path. To visualize this, we plotted trails on the leaflet map on a topographic imagery map layer, a satellite imagery map layer and an open street imagery layer, and added the control to switch between those three imagery layers.

#calculate each part of the length of the trail
length <- 
  macleish_layers[["trails"]] %>%
  group_by(name) %>%
  st_length()

#make a new column in newTrails data frame
newTrails <- 
  newTrails %>%
  mutate(length = 1)

#change the value in the length column to the corresponding value in length set
for (k in 1:nrow(newTrails)){
  newTrails$length[k] <- length[k]
}

#get the total length of each trail
newTrails <- 
  newTrails %>%
  group_by(name) %>% 
  summarise(length = sum(length)) %>%
  arrange(desc(length))
trailInter <- 
  elevation %>%
  st_intersection(macleish_layers[["trails"]])

trailInter <- 
  trailInter %>%
  group_by(name) %>%
  summarise(ele = ifelse(n_distinct(CONTOUR_FT) != 1, sd(CONTOUR_FT), 0)) %>%
  arrange(desc(ele))

We want to look at the exact number of the length and the change in elevation to get a more accurate rating. Therefore, we produced the table newTrails calculate the length of each path and the table trailInter to calculate the change in elevation. For the elevation change, we figured the standard deviation of the elevation of each trail. High standard deviation represents the higher change in elevation and more difficulty, while low standard deviation represents less change in elevation and less difficulty. If a trail doesn’t show up in our table, we regarded it as there is no intersection between the contours lines and that trail which means barely any change in elevation.

Looking at the two tables, we found that Snowmobile Trail and Western Loop are the two largest change in elevation trails and the two longest trails of the all. Therefore, we classified them as in the level of “Difficult.” Path Easy Out, Driveway, and entry trail are the three shortest trails. Since they don’t show up in the change in elevation table, we believed that those trails have the least changes in elevation. Therefore, we rated these three trails as in the level of “Easy.” For the rest of the trails, including Western Loop, Poplar Hill Road, Porcupine Trail, and Vernal Pool Loop, we rated them in the level of “Moderate,”

joinTrail <- 
  macleish_layers[["trails"]] %>%
  st_join(trailInter, join = st_crosses)
## although coordinates are longitude/latitude, st_crosses assumes that they are planar
newTrailsNone <- 
  st_set_geometry(newTrails, NULL)
trailInterNone <- 
  st_set_geometry(trailInter, NULL)

joined <- 
  newTrailsNone %>%
  left_join(trailInterNone, by = "name")


joined$ele[is.na(joined$ele)] <- 0

joined <- 
  joined %>%
  mutate(difficulty = 0.4 * as.numeric(length) + 0.6 * as.numeric(ele)) %>%
  arrange(desc(difficulty))


joined$ID <- seq.int(nrow(joined))

joined <- 
  joined %>%
  mutate(level = ifelse(ID < 4, "hard", ifelse(ID > 6, "easy", "moderate")))
macleish_layers[["trails"]] <- 
  macleish_layers[["trails"]] %>%
  inner_join(joined, by = "name")
pall <- colorNumeric(palette="viridis",
                     domain= macleish_layers[["contours_3m"]]$ELEV_M)
palTrail <- colorFactor(c("red", "blue", "black"), domain = macleish_layers[["trails"]]$level)

leaflet(data=macleish_layers[["trails"]]) %>%
  addTiles() %>%
  addPolylines(data=macleish_layers[["contours_3m"]], 
               color = ~pall(ELEV_M), 
               group = "10'Contours", weight = 1) %>%
  addLegend("bottomleft", pal = pall, 
            values = ~macleish_layers[["contours_3m"]]$ELEV_M, 
            title = "Macleish Contour Line in meters") %>%
  addPolylines(data = macleish_layers[["trails"]], weight = 2, label = ~name, color = ~palTrail(level)) %>%
  addLegend("bottomright", pal = palTrail, 
            values = ~macleish_layers[["trails"]]$level, 
            title = "Difficulty levels for biking") %>%
  addMeasure() %>%
  addControl("Macleish Trails Ranked by Difficulty Levels", position = "topleft") %>%
  addMiniMap(
    toggleDisplay = TRUE) %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Topo") %>% 
  # ...and a satellite imagery map layer
  addProviderTiles("Esri.WorldImagery", group = "Sat") %>%
  # Next you want to add the control to switch between the two
  addLayersControl(baseGroups = c("OpenSreetMap", "Topo","Sat"),
                   options = layersControlOptions(collapsed = F)) %>%
  addEasyButton(easyButton(
    icon="fa-crosshairs", title="Locate Me",
    onClick=JS("function(btn, map){ map.locate({setView: true}); }")))

  1. Github link: https://github.com/ChristinaLyu/miniProject3.git↩